* local projections DID
* lag, single diff
* two lag, single diff
* clean controls may include future treated. 
* reweighted, or other approach? 

clear all
import sas using "/projects/users/########/Snapshot2022/IntermediateData/all_cntiesxmile_windwrkrs2022.sas7bdat", case(lower)
keep if !missing(_id_)
drop if year<2000
drop if year>2020

sort ctyfips year
egen long county = group(ctyfips)
sort county year
xtset county year

* capacity in MW (currently in kW)
gen ccap = c_tcap/1000
gen ccap20 = c_tcap20/1000
gen ccap40 = c_tcap40/1000
gen ccap60 = c_tcap60/1000
gen ccap80 = c_tcap80/1000
gen ccap100 = c_tcap100/1000

bysort county (year): egen mxcap = max(ccap)
sort county year

* outcome variables
gen asin_avg = asinh(avg_emp_earn)
gen asin_white = asinh(white_emp_avgearn)
gen asin_black = asinh(black_emp_avgearn)
gen asin_indnat = asinh(indnat_emp_avgearn)
gen asin_hisp = asinh(hisp_emp_avgearn)
gen asin_male = asinh(male_emp_avgearn)
gen asin_female = asinh(female_emp_avgearn)
gen asin_nohigh = asinh(nohigh_emp_avgearn)
gen asin_highsch = asinh(highsch_emp_avgearn)
gen asin_somecoll = asinh(somecoll_emp_avgearn)
gen asin_college = asinh(college_emp_avgearn)

* consider treated if capacity in county.  
bysort county: egen mxcapc = max(ccap)
bysort county: egen mxcap20 = max(ccap20)
/*
bysort county: egen mxcap40 = max(ccap40)
bysort county: egen mxcap60 = max(ccap60)
bysort county: egen mxcap80 = max(ccap80)
bysort county: egen mxcap100 = max(ccap100)
*/
sort county year

gen grpc = 0
replace grpc = 1 if mxcapc>10 & !missing(mxcapc)
gen grp20 = 0
replace grp20 = 1 if mxcap20>10 & !missing(mxcap20)
/*
gen grp40 = 0
replace grp40 = 1 if mxcap40>10 & !missing(mxcap40)
gen grp60 = 0
replace grp60 = 1 if mxcap60>10 & !missing(mxcap60)
gen grp80 = 0
replace grp80 = 1 if mxcap80>10 & !missing(mxcap80)
gen grp100 = 0
replace grp100 = 1 if mxcap100>10 & !missing(mxcap100)
*/

gen trtc = 0
replace trtc = 1 if grpc==1 & ccap>10 & !missing(ccap)
gen trt20 = 0
replace trt20 = 1 if grp20==1 & ccap20>10 & !missing(ccap20)
/*
gen trt40 = 0
replace trt40 = 1 if grp40==1 & ccap40>10 & !missing(ccap40)
gen trt60 = 0
replace trt60 = 1 if grp60==1 & ccap60>10 & !missing(ccap60)
gen trt80 = 0
replace trt80 = 1 if grp80==1 & ccap80>10 & !missing(ccap80)
gen trt100 = 0
replace trt100 = 1 if grp100==1 & ccap100>10 & !missing(ccap100)
*/

*gen stfips = substr(ctyfips,1,2)
egen state = group(stfips)
egen styr = group(stfips year)

* convert to GW
gen plantcapc = ccap/1000
gen plantcap20 = ccap20/1000
gen plantcap40 = ccap40/1000
gen plantcap60 = ccap60/1000
gen plantcap80 = ccap80/1000
gen plantcap100 = ccap100/1000

gen dcap40 = plantcap40 - plantcap20

gen dcap60 = plantcap60 - plantcap40

gen dcap80 = plantcap80 - plantcap60

gen dcap100 = plantcap100 - plantcap80



* make sure event_date = . if never treated. 

by county (year), sort: gen byte first = trtc==1
by county (year), sort: gen byte first2 = sum(first)==1
gen trtyr = 0 
replace trtyr = year if first2==1
bysort county (year): egen event_date = max(trtyr) 
bysort county (year): replace event_date = . if grpc==0

gen event_time = year - event_date 

global post_window 6
global pre_window 6


xtset county year


/*
Do 4 versions: lag0, with spatial lags (4)
- for discrete and continuous treatment,
- for employment and asin_bestwages if unemp_instate<0.75
*/


* LP-DID Event Study

capture postutil clear
tempfile lpdidc_estimates2022
postfile handle str32 vble double ///
	b_empt_avg vnc_empt_avg b_empc_avg vnc_empc_avg ///
	b_wgt_avg vnc_wgt_avg b_wgc_avg vnc_wgc_avg ///
	b_empt_white vnc_empt_white b_empc_white vnc_empc_white ///
	b_wgt_white vnc_wgt_white b_wgc_white vnc_wgc_white ///
	b_empt_black vnc_empt_black b_empc_black vnc_empc_black ///
	b_wgt_black vnc_wgt_black b_wgc_black vnc_wgc_black ///
	b_empt_indnat vnc_empt_indnat b_empc_indnat vnc_empc_indnat ///
	b_wgt_indnat vnc_wgt_indnat b_wgc_indnat vnc_wgc_indnat ///
	b_empt_hisp vnc_empt_hisp b_empc_hisp vnc_empc_hisp ///
	b_wgt_hisp vnc_wgt_hisp b_wgc_hisp vnc_wgc_hisp ///
	b_empt_male vnc_empt_male b_empc_male vnc_empc_male ///
	b_wgt_male vnc_wgt_male b_wgc_male vnc_wgc_male ///
	b_empt_female vnc_empt_female b_empc_female vnc_empc_female ///
	b_wgt_female vnc_wgt_female b_wgc_female vnc_wgc_female ///
	b_empt_nohigh vnc_empt_nohigh b_empc_nohigh vnc_empc_nohigh ///
	b_wgt_nohigh vnc_wgt_nohigh b_wgc_nohigh vnc_wgc_nohigh ///
	b_empt_highsch vnc_empt_highsch b_empc_highsch vnc_empc_highsch ///
	b_wgt_highsch vnc_wgt_highsch b_wgc_highsch vnc_wgc_highsch ///
	b_empt_somecoll vnc_empt_somecoll b_empc_somecoll vnc_empc_somecoll ///
	b_wgt_somecoll vnc_wgt_somecoll b_wgc_somecoll vnc_wgc_somecoll ///
	b_empt_college vnc_empt_college b_empc_college vnc_empc_college ///
	b_wgt_college vnc_wgt_college b_wgc_college vnc_wgc_college ///
	using `lpdidc_estimates2022', replace


foreach x in avg white black indnat hisp male female nohigh highsch somecoll college {

	* forward changes in outcomes, two years prior as base year: 
	* gen forward changes - employ
	forval j = 0/$post_window {
		gen D`j'y = F`j'.`x'_emp_in - L2.`x'_emp_in
	}
	gen Dn1y = L.`x'_emp_in - L2.`x'_emp_in
	forval j = 3/$pre_window {
		gen Dm`j'y = L`j'.`x'_emp_in - L2.`x'_emp_in
	}
	* gen forward changes - log earnings
	forval j = 0/$post_window {
		gen D`j'ya = F`j'.asin_`x' - L2.asin_`x'
	}
	gen Dn1ya = L.asin_`x'  - L2.asin_`x' 
	forval j = 3/$pre_window {
		gen Dm`j'ya = L`j'.asin_`x'  - L2.asin_`x' 
	}

	* 1. dummy treat, single diff treat, no lag, not treated as of t+h, not reweighted, surrounding capacity controls.
	* employ_in
	matrix blpdid = J(1,13,.)
	matrix blpdidvar = J(1,13,.)
	matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	forval j = 0/$post_window {
		reghdfe D`j'y D.trtc D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
			if D.trtc==1 | F`j'.trtc==0,  ///
			absorb(styr) vce(cluster county)
		matrix blpdid[1,`j'+$post_window + 1] = _b[D.trtc]
		matrix blpdidvar[1,`j'+$post_window +1] = sqrt(e(V)[1,1])
		
		if `j'>2 & `j'<=$pre_window {
			reghdfe Dm`j'y D.trtc D.dcap40 D.dcap60 D.dcap80 D.dcap100 ///
				if D.trtc==1 | trtc==0, ///
				absorb(styr) vce(cluster county)
			local p = $post_window +1 -`j'
			matrix blpdid[1,`p'] = _b[D.trtc]
			matrix blpdidvar[1,`p'] = sqrt(e(V)[1,1])
		}
	} 
	reghdfe Dn1y 			///
		D.trtc D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// 
		if D.trtc == 1 | trtc ==0, 	///  /* clean controls */
		absorb(styr) vce(cluster county)	
	matrix blpdid[1,$post_window ] = _b[D.trtc]
	matrix blpdidvar[1,$post_window ] = sqrt(e(V)[1,1])

	matrix blpdid[1,$post_window -1] = 0
	matrix blpdidvar[1,$post_window -1] = 0	

	matrix blp_empt_`x' = (blpdid \ blpdidvar)
	matrix rownames blp_empt_`x' = b_empt_`x' var_empt_`x'
	matrix drop blpdid 
	matrix drop blpdidvar

	matrix B_`x' = (blp_empt_`x')

	* asin_bestwages
	matrix blpdid = J(1,13,.)
	matrix blpdidvar = J(1,13,.)
	matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	forval j = 0/$post_window {
		reghdfe D`j'ya D.trtc D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
			if D.trtc==1 | F`j'.trtc==0,  ///
			absorb(styr) vce(cluster county)
		matrix blpdid[1,`j'+$post_window + 1] = _b[D.trtc]
		matrix blpdidvar[1,`j'+$post_window +1] = sqrt(e(V)[1,1])
		
		if `j'>2 & `j'<=$pre_window {
			reghdfe Dm`j'ya D.trtc D.dcap40 D.dcap60 D.dcap80 D.dcap100	 ///
				if D.trtc==1 | trtc==0, ///
				absorb(styr) vce(cluster county)
			local p = $post_window +1 -`j'
			matrix blpdid[1,`p'] = _b[D.trtc]
			matrix blpdidvar[1,`p'] = sqrt(e(V)[1,1])
		}
	} 
	reghdfe Dn1ya 			///
		D.trtc  D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// /* lags */
		if D.trtc == 1 | trtc ==0, 	///  /* clean controls */
		absorb(styr) vce(cluster county)	
	matrix blpdid[1,$post_window ] = _b[D.trtc]
	matrix blpdidvar[1,$post_window ] = sqrt(e(V)[1,1])

	matrix blpdid[1,$post_window -1] = 0
	matrix blpdidvar[1,$post_window -1] = 0	

	matrix blp_wgt_`x' = (blpdid \ blpdidvar)
	matrix rownames blp_wgt_`x'  = b_wgt_`x' var_wgt_`x'
	matrix drop blpdid 
	matrix drop blpdidvar

	matrix B_`x' = (B_`x' \ blp_wgt_`x')


	************************
	* Continuous treatment *
	************************

	* employ_in
	matrix blpdid = J(1,13,.)
	matrix blpdidvar = J(1,13,.)
	matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	forval j = 0/$post_window {
		reghdfe D`j'y D.plantcapc D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
			if D.plantcapc>=0.01 | F`j'.plantcapc==0,  ///
			absorb(styr) vce(cluster county)
		matrix blpdid[1,`j'+$post_window + 1] = _b[D.plantcapc]
		matrix blpdidvar[1,`j'+$post_window +1] = sqrt(e(V)[1,1])
		
		if `j'>2 & `j'<=$pre_window {
			reghdfe Dm`j'y D.plantcapc D.dcap40 D.dcap60 D.dcap80 D.dcap100 ///
				if D.plantcapc>=0.01 | plantcapc==0, ///
				absorb(styr) vce(cluster county)
			local p = $post_window +1 -`j'
			matrix blpdid[1,`p'] = _b[D.plantcapc]
			matrix blpdidvar[1,`p'] = sqrt(e(V)[1,1])
		}
	} 
	reghdfe Dn1y 			///
		D.plantcapc D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// 
		if D.plantcapc >=0.01 | plantcapc ==0, 	///  /* clean controls */
		absorb(styr) vce(cluster county)	
	matrix blpdid[1,$post_window ] = _b[D.plantcapc]
	matrix blpdidvar[1,$post_window ] = sqrt(e(V)[1,1])

	matrix blpdid[1,$post_window -1] = 0
	matrix blpdidvar[1,$post_window -1] = 0	

	matrix blp_empc_`x' = (blpdid \ blpdidvar)
	matrix rownames blp_empc_`x'  = b_empc_`x' var_empc_`x'
	matrix drop blpdid 
	matrix drop blpdidvar

	matrix B_`x' = (B_`x' \ blp_empc_`x')

	* asin_bestwages
	matrix blpdid = J(1,13,.)
	matrix blpdidvar = J(1,13,.)
	matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
	forval j = 0/$post_window {
		reghdfe D`j'ya D.plantcapc 	D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
			if D.plantcapc>=0.01 | F`j'.plantcapc==0,  ///
			absorb(styr) vce(cluster county)
		matrix blpdid[1,`j'+$post_window + 1] = _b[D.plantcapc]
		matrix blpdidvar[1,`j'+$post_window +1] = sqrt(e(V)[1,1])
		
		if `j'>2 & `j'<=$pre_window {
			reghdfe Dm`j'ya D.plantcapc D.dcap40 D.dcap60 D.dcap80 D.dcap100	 ///
				if D.plantcapc>=0.01 | plantcapc==0, ///
				absorb(styr) vce(cluster county)
			local p = $post_window +1 -`j'
			matrix blpdid[1,`p'] = _b[D.plantcapc]
			matrix blpdidvar[1,`p'] = sqrt(e(V)[1,1])
		}
	} 
	reghdfe Dn1ya 			///
		D.plantcapc  D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// /* lags */
		if D.plantcapc >=0.01 | plantcapc ==0, 	///  /* clean controls */
		absorb(styr) vce(cluster county)	
	matrix blpdid[1,$post_window ] = _b[D.plantcapc]
	matrix blpdidvar[1,$post_window ] = sqrt(e(V)[1,1])

	matrix blpdid[1,$post_window -1] = 0
	matrix blpdidvar[1,$post_window -1] = 0	

	matrix blp_wgc_`x' = (blpdid \ blpdidvar)
	matrix rownames blp_wgc_`x'  = b_wgc_`x' var_wgc_`x'
	matrix drop blpdid 
	matrix drop blpdidvar

	matrix B_`x' = (B_`x' \ blp_wgc_`x')

	drop D0y-Dm6ya

}



foreach v in pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6 {
	post handle ("`v'") ///
	(B_avg["b_empt_avg", "`v'"]) (B_avg["var_empt_avg", "`v'"]) ///
	(B_avg["b_empc_avg", "`v'"]) (B_avg["var_empc_avg", "`v'"]) ///
	(B_avg["b_wgt_avg", "`v'"]) (B_avg["var_wgt_avg", "`v'"]) ///
	(B_avg["b_wgc_avg", "`v'"]) (B_avg["var_wgc_avg", "`v'"]) ///
	(B_white["b_empt_white", "`v'"]) (B_white["var_empt_white", "`v'"]) ///
	(B_white["b_empc_white", "`v'"]) (B_white["var_empc_white", "`v'"]) ///
	(B_white["b_wgt_white", "`v'"]) (B_white["var_wgt_white", "`v'"]) ///
	(B_white["b_wgc_white", "`v'"]) (B_white["var_wgc_white", "`v'"]) ///
	(B_black["b_empt_black", "`v'"]) (B_black["var_empt_black", "`v'"]) ///
	(B_black["b_empc_black", "`v'"]) (B_black["var_empc_black", "`v'"]) ///
	(B_black["b_wgt_black", "`v'"]) (B_black["var_wgt_black", "`v'"]) ///
	(B_black["b_wgc_black", "`v'"]) (B_black["var_wgc_black", "`v'"]) ///
	(B_indnat["b_empt_indnat", "`v'"]) (B_indnat["var_empt_indnat", "`v'"]) ///
	(B_indnat["b_empc_indnat", "`v'"]) (B_indnat["var_empc_indnat", "`v'"]) ///
	(B_indnat["b_wgt_indnat", "`v'"]) (B_indnat["var_wgt_indnat", "`v'"]) ///
	(B_indnat["b_wgc_indnat", "`v'"]) (B_indnat["var_wgc_indnat", "`v'"]) ///
	(B_hisp["b_empt_hisp", "`v'"]) (B_hisp["var_empt_hisp", "`v'"]) ///
	(B_hisp["b_empc_hisp", "`v'"]) (B_hisp["var_empc_hisp", "`v'"]) ///
	(B_hisp["b_wgt_hisp", "`v'"]) (B_hisp["var_wgt_hisp", "`v'"]) ///
	(B_hisp["b_wgc_hisp", "`v'"]) (B_hisp["var_wgc_hisp", "`v'"]) ///
	(B_male["b_empt_male", "`v'"]) (B_male["var_empt_male", "`v'"]) ///
	(B_male["b_empc_male", "`v'"]) (B_male["var_empc_male", "`v'"]) ///
	(B_male["b_wgt_male", "`v'"]) (B_male["var_wgt_male", "`v'"]) ///
	(B_male["b_wgc_male", "`v'"]) (B_male["var_wgc_male", "`v'"]) ///
	(B_female["b_empt_female", "`v'"]) (B_female["var_empt_female", "`v'"]) ///
	(B_female["b_empc_female", "`v'"]) (B_female["var_empc_female", "`v'"]) ///
	(B_female["b_wgt_female", "`v'"]) (B_female["var_wgt_female", "`v'"]) ///
	(B_female["b_wgc_female", "`v'"]) (B_female["var_wgc_female", "`v'"]) ///
	(B_nohigh["b_empt_nohigh", "`v'"]) (B_nohigh["var_empt_nohigh", "`v'"]) ///
	(B_nohigh["b_empc_nohigh", "`v'"]) (B_nohigh["var_empc_nohigh", "`v'"]) ///
	(B_nohigh["b_wgt_nohigh", "`v'"]) (B_nohigh["var_wgt_nohigh", "`v'"]) ///
	(B_nohigh["b_wgc_nohigh", "`v'"]) (B_nohigh["var_wgc_nohigh", "`v'"]) ///
	(B_highsch["b_empt_highsch", "`v'"]) (B_highsch["var_empt_highsch", "`v'"]) ///
	(B_highsch["b_empc_highsch", "`v'"]) (B_highsch["var_empc_highsch", "`v'"]) ///
	(B_highsch["b_wgt_highsch", "`v'"]) (B_highsch["var_wgt_highsch", "`v'"]) ///
	(B_highsch["b_wgc_highsch", "`v'"]) (B_highsch["var_wgc_highsch", "`v'"]) ///
	(B_somecoll["b_empt_somecoll", "`v'"]) (B_somecoll["var_empt_somecoll", "`v'"]) ///
	(B_somecoll["b_empc_somecoll", "`v'"]) (B_somecoll["var_empc_somecoll", "`v'"]) ///
	(B_somecoll["b_wgt_somecoll", "`v'"]) (B_somecoll["var_wgt_somecoll", "`v'"]) ///
	(B_somecoll["b_wgc_somecoll", "`v'"]) (B_somecoll["var_wgc_somecoll", "`v'"]) ///
	(B_college["b_empt_college", "`v'"]) (B_college["var_empt_college", "`v'"]) ///
	(B_college["b_empc_college", "`v'"]) (B_college["var_empc_college", "`v'"]) ///
	(B_college["b_wgt_college", "`v'"]) (B_college["var_wgt_college", "`v'"]) ///
	(B_college["b_wgc_college", "`v'"]) (B_college["var_wgc_college", "`v'"]) 
}

postclose handle
preserve
use "/projects/users/########/Snapshot2022/Results/lpdid2_strateq2_county2.dta", clear
append using `lpdidc_estimates2022'
save "/projects/users/########/Snapshot2022/Results/lpdid2_strateq2_county2.dta", replace
restore

